#include "Analyze.h"


Analyze::Analyze(string fl1)
{
    REF= new Reader(fl1);
    LD = REF->LibelaDat();
    REF->Logreader();
    T = REF->temp;
}



vector<long double> Analyze::MinFinder(){


    long double min= LD[0][1];
    long double ic = 0;
    int nlines = LD.size();
    for (int k1=0; k1<nlines; k1++){
        if (LD[k1][1]<=min){
            min = LD[k1][1];
            ic = k1;
        }

    }
    vector<long double> fin;
    fin.push_back(ic);
    fin.push_back(min);
    return fin;
}

vector<long double> Analyze::Correlation(int optdim, int trunk){

        vector<long double> corr;
        int k, t;
        int corrsize = LD.size()-1;
        double tck, norm;
        for (t=0; t<corrsize; t++){
            norm = norm + ((LD[t][optdim])*(LD[t][optdim]))/(LD.size()-1);
        }
        for (k=0; k<(trunk); k++){
            tck=0;
            for (t=0; t<(corrsize-k); t++){
                tck = tck + ((LD[t][optdim])*(LD[t+k][optdim]))/((LD.size()-1-k)*norm);
            }

            corr.push_back(abs(tck));

        }
        return corr;
}

vector<long double> Analyze::BoltzmannWeightConvergence(int trunk){
    long double norm;
    long double tener;
    vector< vector<long double> > BOLT;
    vector <long double > tbolt;
    int size = LD.size();
    if (size<trunk){
        cout << "Trunk value > Simulation Steps"<< " (Steps = " << LD.size() << ")" << endl;
    }
    else{

    for (int k1=0; k1<trunk; k1++){
                tener=0;
                norm =0;
                for (int k2=0; k2<=k1; k2++){
                            norm = norm + exp(-LD[k2][1]/(K*T));
                            }
                for (int k2=0; k2<=k1; k2++){
                            tener = tener + (LD[k2][1])*exp(-LD[k2][1]/(K*T))/norm;

                    }

                tbolt.push_back(tener);

    }
    return tbolt;
    }
}

// Calculate Internal energy of the system
long double Analyze::MeanInternalEnergy(){
    long double norm;
    long double tener;
    int size = LD.size();
    for (int k1=0; k1<size; k1++){
                tener=0;
                norm =0;
                for (int k2=0; k2<=k1; k2++){
                            norm = norm + exp(-LD[k2][1]/(K*T));
                            }
                for (int k2=0; k2<=k1; k2++){
                            tener = tener + (LD[k2][1])*exp(-LD[k2][1]/(K*T))/norm;

                    }
             }
    cout << tener << endl;
    return tener;
    }

//Calculates the Mean Potential non-covalent energy
long double Analyze::Mean(){

    long double tmedian = 0;
    int nlines = LD.size();
    for (int k1=0; k1<nlines; k1++){
            tmedian = tmedian + (LD[k1][1])/(nlines);
        }

    return tmedian;
}




//Calculates the convergence fo the mean.
 vector <long double >  Analyze::MeanConvergence(){
        long double tmedian;
        vector <long double > Meanvec;
        vector< vector <long double > > Result;
        int nlines = LD.size();
        for (int k1=0; k1<nlines; k1++){
            tmedian=0;
            for (int k2=0; k2<=k1; k2++){
                tmedian = tmedian + (LD[k2][1])/(k1+1);
            }
            Meanvec.push_back(tmedian);
        }
        return Meanvec;
}


// normaliza e a anliza pelo minimo
vector< long double > Analyze::NormByMin(){

            vector <long double > Normali;
            long double texp;
            long double min = this->MinFinder()[1];
            int nlines = LD.size();
            for (int k1=0; k1<nlines; k1++){
                texp = -(LD[k1][1]-min);
                Normali.push_back(texp);

            }
            return Normali;
            }





//Free Energy Calculator
long double Analyze::FreeEnergy(){
        long double tener = 0;
        int nlinesEnergy = LD.size();


                    long double norm =0;
                    long double Angular_Factor = 0;
                    long double Translational_Factor = 0;
                    long double FE;

                    int rNconf = REF->Rnconf;

                    for (int i=0; i<rNconf-1; i++){
                            norm = norm + exp(-((REF->enerconf[i])/(K*T)));

                        }

                    for (int i=0; i<nlinesEnergy; i++){
                                Angular_Factor = (abs(LD[i][6])*abs(LD[i][7])*abs(LD[i][8])*PI*PI*PI)/(8*PI*PI*180*180*180);
                                Translational_Factor = (abs(LD[i][3])*abs(LD[i][4])*abs(LD[i][5]))/(REF->volx*REF->voly*REF->volz);
                                tener = tener + exp(-(LD[i][1]+LD[i][10]/(K*T)))*Angular_Factor*Translational_Factor/(norm);

                    }
                    FE = -K*T*log(tener);
                    cout << FE << endl;
                    return FE;
}

vector <long double> Analyze::ConvFreeEnergy(){

        long double tener = 0;
        int nlinesEnergy = LD.size();


                    long double norm =0;
                    long double Angular_Factor = 0;
                    long double Translational_Factor = 0;
                    vector < long double> conv;
                    long double FE;
                    int rNconf = REF->Rnconf;
                    for (int i=0; i<=rNconf; i++){

                            norm = norm + exp(-((REF->enerconf[i])/(K*T)));


                        }

                    for (int i=0; i<nlinesEnergy; i++){
                                Angular_Factor = (abs(LD[i][6])*abs(LD[i][7])*abs(LD[i][8])*PI*PI*PI)/(8*PI*PI*180*180*180);
                                Translational_Factor = (abs(LD[i][3])*abs(LD[i][4])*abs(LD[i][5]))/(REF->volx*REF->voly*REF->volz);
                                tener = tener + exp(-(LD[i][1]+LD[i][10]/(K*T)))*Angular_Factor*Translational_Factor/(norm);
                                FE = -K*T*log(tener);
                                conv.push_back(FE);

                    }
                    return conv;
}


long double Analyze::Entropy(){
    long double Ent = -(this->FreeEnergy()-this->MeanInternalEnergy())/this->T;
    return Ent;
}


